using QuantumOptics, OrdinaryDiffEq, BenchmarkTools

b = SpinBasis(1//2)
psi0 = spindown(b)
u0 = ComplexF64[0.5, 0.75] 
sc = semiclassical.State(psi0, u0)
t₀, t₁ = (0.0, pi)
const σx = sigmax(b)

fquantum(t, q, u) = σx + cos(u[1])*identityoperator(σx)
fclassical!(du, u, q, t) = (du[1] = sin(u[2]); du[2] = 2*u[1])
f!(dstate, state, p, t) = semiclassical.dschroedinger_dynamic!(dstate, fquantum, fclassical!, state, t)
prob = ODEProblem(f!, sc, (t₀, t₁))

julia> @benchmark sol = solve($prob, DP5(); save_everystep=false)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  50.166 μs … 133.706 ms  ┊ GC (min … max):  0.00% … 99.93%
 Time  (median):     52.583 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   70.621 μs ±   1.339 ms  ┊ GC (mean ± σ):  23.35% ±  3.73%

       ▃█▅                                                      
  ▁▂▃▆████▇▃▂▂▂▂▂▂▂▂▂▃▃▃▃▂▃▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  50.2 μs         Histogram: frequency by time         69.8 μs <

 Memory estimate: 93.81 KiB, allocs estimate: 1136.

julia> @benchmark semiclassical.schroedinger_dynamic([$t₀, $t₁], $sc, $fquantum, $fclassical!)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  117.875 μs … 137.310 ms  ┊ GC (min … max):  0.00% … 99.87%
 Time  (median):     122.959 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   148.349 μs ±   1.379 ms  ┊ GC (mean ± σ):  15.47% ±  6.25%

    ▃█▇▁▁▃▃                                                      
  ▂▅████████▄▃▃▂▂▃▃▂▃▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  118 μs           Histogram: frequency by time          163 μs <

 Memory estimate: 225.66 KiB, allocs estimate: 2637.